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SUMMARY 

The study of a thin, incompressible Newtonian fluid layer trapped between two 
almost parallel, sliding surfaces has been actively pursued in the last decades. This 
subject includes lubrication applications such as slider bearings or the sealing of non- 
pressurized fluids with rubber rbtary shaft seals. In the present work we analyze 
numerically the flow of lubricant fluid through a micro-gap of sealing devices. The first 
stage of this study is carried out assuming that a ‘small-gap’ parameter 6 attains an 
extreme value in the Navier-Stokes equations. The precise meaning of small-gap is 
achieved by the particular limit 6 = 0 which, within the bounds of the hypotheses, 
predicts transport of lubricant through the sealed area by centrifugal instabilities. 
Numerical results obtained with the penalty function approximation in the finite 
element method are presented. In particular, the influence of inflow and outflow 
boundary conditions, and their impact in the simulated flow are discussed. 

INTRODUCTION 

Most seals are relatively simple elements widely employed in diverse types of 
rotary machines. This sealing component is often used to seal rotating shafts at low oil 
pressures, avoiding the transport of cont ami nant to, or lubricant from, the equipment it 
protects. The seal, bonded to the oil reservoir, is stationary and presents a narrow 
section that slides over the moving surface of the rotary shaft. Fig.l shows a cross 
section of the geometry under consideration. 

The device is designed to have an interference with the shaft. Therefore, once the 
piece is mounted, the compliance of the elaistic body ensures a perfect fit between the 
seail and the cylindricad surface of the shaft. Under these conditions, some of the initiad 
asperities of the sead wear out after a brief period of time, leaving an extremely thin 
layer of lubricant fluid that separates the surfaces in contact. This was first noticed by 
Jagger (ref.l) and, ever since, numerous explanations attempted to account for two 
consequences of this experimented fact: the hydrodynamic force able to sustain a gap 
between the two bodies and the mechanisms that prevent the fluid from leaking 
through. Jagger proposed that the surface tension of the seeded fluid controls leakage 
thanks to a meniscus formed on the air side. Kawahara and Hirabayashi (ref.2) observed 
that a properly installed and functionad seal leaked when the installation is reversed. 

With the assumption of relative paradlel sliding between two rough surfaces, 
lubrication theory hais been the chosen tool by many researchers to answer these 
fundamental problems (see e.g. ref.3). The load-cauxying capacity of parallel sliding of 
rough surfaces was first studied by Davies (ref.4). Later on, Jagger and Walker (ref.5) 
aissumed that the asperities act as micro-bearings pads in the contact area. However, 
Lebeck (refs. 6 and 7) concluded that none of the existing models can fully explain the 
sliding motion as commonly observed in experiments. Gabelli and Poll (ref. 8) studied 
the dominant action of the surfatce microgeometry in the formation of the lubricant 
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film. They found that the contribution of mechanical pressure to the load- carrying 
capacity due to body contact is very small and indeed negligible. Sal ant (ref.9) claimed 
that micro-undulations in the lip surface restrict leakage by virtue of a ‘reverse- 
pumping’ process in which fluid is driven from the low to the high pressure side. 
However, no one has really observed such micro- undulations, either in static or dynamic 
conditions (ref. 10). 

Combinations of angular velocity and system eccentricity beyond the ability of 
the sealing device to maintain contact with the shaft would cause the seal to leak 
profusely. It has been suggested that an inherent pumping mechanism (refill), sufficient 
to counterbalance those influences promoting leakage, would be given by a relative 
motion between the sealing surfaces. Besides all these explanations, at present there is a 
wide gap between theory and practice, and a feasible explanation of the mechanism 
involved in the seeding action is still pending, even though elastomeric seals have been 
used extensively since the 1940’s. _ . , . 

In the next section, we establish the small-gap equations using a rather simple 
order-of-magnitude analysis. This is followed by numerical examples showing the 
validity of the proposed model and the influence of the boundary conditions in the 
numerical predictions. 



ANALYTICAL MODEL 

We assume an oil-film already formed ignoring any mechanical contact between 
the sealing device and the shaft, as well as any distortion of the upper elastic seal. We 
consider a thin viscous liquid layer bounded above by a smooth surface and below by a 
perfectly rounded shaft, without including edge effects such as the meniscus 
experimentally observed on the air-side (ref.l). Despite the fact that the film within the 
gap is very thin, we assume it to be thick enough to conform to the continuum 
hypothesis. There is no local rupture of the film such as cavitation or dry spots in the 
contact area, and the layer consists of an incompressible Newtonian fluid with constant 
properties under isothermal conditions. 

We begin with the Navier-Stokes equations written in cylindrical coordinates 
(ref.12), setting the direction of the line r = 0 coincident with the shaft axis 
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In absence of a free surface the gravitational body force is expressed as the 
gradient of a scalar quantity and, therefore, it has been included in the pressure 
gradient term. 

The analysis of the lubricant flow in the micro-gap involves, roughly speaking, 
three disparate length scales, namely, the radius R of the shaft ( ~ 0.04m), the much 
smaller thickness h 0 of the fluid ( ~ 10 pm) and an intermediate length b characterizing 
the axial extent of the contact region ( ~ 200 pm) (see fig. 1). Next, the Navier-Stokes 
equations, once recasted with the aforementioned scales, are simplified by letting the 
ratio between the gap height and the radius of the shaft formally approach zero. 


Inner region: lubrication regime 


Given the tiny thickness of the lubricant film, radial oscillations proportional to 
the gap height will alter considerably the flow inside the micro-gap. To analyze this 
effect, consider a shaft rotating at angular velocity Cl and separated an average distance 
h a from the stationary seal (fig.l). Introducing now 
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into the equations of motion (l)-(4), and letting $ =-& 
holding everything else fixed, we get ** 


formally approach zero while 
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and 
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R e = — jr 2 : Reynolds number 
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Severed other scaling are possible (ref.13), but this particular choice seems to be 
consistent with Gabelli and Poll observations (ref.8). They stated that the average 
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pressure gradient in the circumferential direction is indeed negligible when compared 
with the pressure gradient across the sealing contact. The squeezing Reynolds number <7 
is commonly so small that inertia terms can be neglected and the classical lubrication 
theory cam be applied. Moreover, for small R,, as it turns out to be in most 
applications, the flow is stable to small perturbations (ref.14). In the absence of 
mechanical vibrations, no secondary flow is possible at this level, the circumferential 
flow being stable and mostly of Couette type. 


Outer region: centrifugal effects 


The loss of contact between seal and shaft, combined with changes in the 
geometry, will introduce different features in the flow behavior as we go farther away 
from the gap. For a slowly-varying channel of height d(z) on the outside (fig.l), we 
rescale the flow field by writing 
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where d is some mean value of d(z). It can be shown that the equations of motion, in 
the limit S = d/R -* 0, become 
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where 
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Note that the above system of equations are the so-called ‘small-gap’ equations, widely 
used in the context of the stability of axisymmetric Taylor-Couette flows (ref.15). 
Again, while the curvature effects are almost completely neglected, they are retained 
through the centrifugal term by imposing the Taylor number be held fixed as 6 -► 0. It 
follows that a rigid seal separated from a rotary shaft by a thin lubricating film is 
subject to centrifugal instabilities in a neighborhood of the contact area, driving 
eventually a secondary flow across the gap. Had we used these scales in the inner region 

— — — " 2 = ( t? J R e 6 -* 0 as<M) (18) 


T a = 


outer region | inner region 


and the Taylor number indicates where curvature effects must be retained, regardless of 
the scales chosen. 

In principle, the flow of lubricant fluid is governed by a set of equations similar 
to those of a stratified flow, where the centripetal acceleration plays the role of the 
buoyancy force, although for rotating flows whose inner surface moves the basic state 
could be unstable to small disturbances. 

The domain under consideration is typically unbounded, and clearly some 
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difficulties arise with the introduction of artificial boundaries. It is known that 
isothermal flows with open boundaries can be successfully modeled employing the so- 
called natural boundary conditions, but other techniques are needed to solve the more 
complicated problem of buoyancy-driven flows, where the use of the natural boundary 
conditions is precluded by the additional pressure gradient generated by the buoyancy 
term (refs. 16 to 18). 

In the next section we address the appropriate use of the penalty formulation in 
the finite element method for unbounded flows in presence of variable body forces. 
Finally, we will see how different open boundary conditions can lead to contradictory 
predictions in the flow behavior. 


PENALTY FUNCTION FORMULATION FOR THE N-S EQUATIONS 

In what follows, we denote the coordinate directions as (x,y) or (r^Xj), the 
transverse velocity components as (u,t>) or (uj,u 2 ), the azimuthal component of the 
velocity as u>, and the pressure as p; 8 { j is the Kronecker delta. For convenience in the 
treatment of the boundary conditions, we rewrite the equations of motion as 
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The Taylor number T a is defined in equation (17), and the stress in the fluid is given 
now by 




( 22 ) 


A weak form is obtained by taking the ijjmer product of the transverse 
momentum equations (20) with a weighting function w = {W v WO, and multiplying the 
azimuthal momentum equation (21) by a scalar function W. The penalty method is 
implemented by introducing the pseudo-constitutive relation (ref. 19) 


P = P ‘~ X W j (23) 

where p t is the hydrostatic pressure for a fluid at rest and A is the penalty parameter. 
Upon application of Green’s theorem and substituting p by the above expression, we get 


fi 3 a 3 3 an 3 

where the surface forces S j and the volume forces V j are defined by 
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(27) 


V,= \t. w 2 6, , W, da 

__ fl 

and 7 * = (n 1 ,n 2 ) is the unit vector normal to the boundary dft and pointing outwards. 
On a vertical boundary with normal n = (1,0), the integrand of Sf reduces to 

W = (W„oy. -p + 2|2 (28) 

for normal traction, and 

w = (o,w,y. + (29) 

for shear traction, in a weak rather than a pointwise sense. 


Boundary conditions 

The boundary conditions are the usual no-slip and no mass penetration at solid 
walls on the physical boundaries. This is, it = v = 0 and w = 1 at the lower boundary 
y = 0, which represents the outer surface of the rotating shaft, and u = v = u; = 0at the 
upper boundary, which is stationary. 

At the open boundaries, on both sides of the contact region, the following two 
open boundary conditions (OBC) are employed for the flow field. 

(i) Stress-free or natural boundary condition (NBC). We set the normal and shear, 
stresses in equation (26) equal to zero. 

(u) Free-boundary condition (FBC). We evaluate the line integrals (26) of the weak 
form of the momentum equations using values computed on the outflow elements. 
Then, we force the line integrals into the right-hand-side of the discretized equations 
until convergence is achieved (refs. 17 and 18). 

The natural boundary condition dw/dn=0 is used in the weak form of the transport 
equation (25) instead. 


FINITE ELEMENT METHOD 

We discretize the domain into M elements and N nodes, and we expand the 
velocity components using bilinear quadrilateral elements and piecewise constant 
elements for pressure. All terms of the weak form of the governing equations are 
evaluated with full Gaussian quadrature, except the penalty term, where selective 
reduced integration is used (ref.20). The weighting functions are set equal to the basis 
functions, except in the convective terms, where perturbed Petrov- Gailer kin functions 
with balancing tensor diffusivity are employed (refs. 21 and 22). The time integration is 
based on the theta method with lumped mass matrices in the time derivatives. The 
numericad evaluation of the weighted residuals of the momentum equations leads to a 
nonlinear system of equations that is solved by Newton iteration using a direct solver 
based on Gauss elimination for unsymmetric banded matrices (ref. 23). A convergence 
tolerance less than 1% of the relative change j| Au 1 ' ||/|| u v || in the velocity field is imposed 
to terminate each full I'-th Newton iteration. The pressure p e over each element fi e is 
calculated using the weak form of the relation (23) 

P e = | V-u dn (30) 
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where the cross bax denotes selective reduced integration. 

To march in time we use the velocity field u n and pressure p n at time t n to 
evaluate the terms of S f and V f of the buoyancy force vector b n . Having determined b n , 
we compute the velocity field using the Newton linearization algorithm. Once u n + 1 is 
known, we update the pressure by means of the equation (30) and solve the transport 
equation for to n + 1 . The scheme is repeated until steady state is achieved. Time 
integration is terminated when the relative change between time steps is 


where 


u n + 1 _ u n 
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for some prescribed error tolerances and e 2 . All the following results are obtained with 
the fully implicit algorithm starting from zero initial conditions. 


NUMERICAL EXAMPLES 

Preliminary computations showed the necessity of using mesh grading as the 
contact region is approached from both sides. Transition elements are also employed to 
avoid extremely small elements in the contact area. 

The geometry and the finite element mesh employed for the present calculations 
axe shown in fig. 2. Details of the discretization can be appreciated in fig. 3. The mesh 
contains 2035 nodes and 1864 elements, and the penalty parameter A is equal to 10 8 in 
all cases. The relevant lengths axe b — 200 pm, which is used as reference length, 
h a = 10 pm, and R = 0.035m. 

The pressure is adjusted at every time step in such a way that is always zero at 
the first element (located at x = — 11, y = 0), and the line integrals (26) axe evaluated, 
if the OBC requires so, with values computed on the outflow elements. 

Results of the transverse velocity field, the azimuthal component of the velocity, 
and the pressure obtained with the FBC at the open boundaries are till shown in fig. 4, 
and continues up to fig. 7. The simulation corresponds to a Taylor number T a = 15. The 
steep pressure gradient developed across the gap is shown in fig.5, and the resultant flux 
of lubricant flowing from the air-side to the oil-side is observed in fig.6. Streamline 
contours axe plotted in fig. 7. Similar results obtained with the NBC at the outlets can 
be seen from fig.8 to fig. 11. The striking differences in the numerical predictions of both 
OBC are better illustrated in fig.7 and fig.ll. The former clearly shows that an 
improper treatment of the open boundary conditions causes backflow into the 
computational domain. The intensity of the returning flow due to the NBC at the air- 
side outflow boundary induces a cell structure in an otherwise alm ost plane Couette 
flow (see w in fig.4). 


CONCLUSIONS 

The geometry and, in particular, the tiny size of the gap imposes a severe 
constraint in the numerical simulation. Furthermore, we have seen that an improper 
specification of the outflow boundary condition can cause artificial returning flow which, 
for the present problem, spoils the solution in the whole computational domain. Both 
boundary conditions show that centrifugal instabilities play an important role in the 
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figure 4 
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figure 8 
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figure 11 
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rnS^,V+ eCtaiU wS^ 0f devices bere studied, even when their predictions are 

achlvi 1^- ^ le ^l NBC u SU . ggeS ^ S leaka S e ’ tbe FBC indicates that sealing^' 

Sd* wh d by PT-r 2 , 0 ? 1 fr °? 1 1 . the air ‘ side > wb ere t^e azimuthal flow is stable, to the^oil- 
side, where centnfugal instabilities set in 

mTrimofthe cUptSnlTdomS (ref^t unwe,come lrfkcti< “ s ‘°™ds 

. Brides its simplicity, the capability of the small-gap limit in incorporating the 

established Other 0 g2 ub , ncant *J uid through the micro-gap of sealing devices has been 
established. Other effects such as capillary forces on the oil-air interface and 
temperature variations should be included in future works. 
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